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O/VP^^^V  i/O.  ABSTRACT 

An  investigation  was  maJ|    fcflfcmckling  of  rectangular  plates 
with  free,  simply  supported,  or  clamped  edges  and  loaded  in  uniform 
compression  or  pure  shear.   The  finite  element  method  with  a  rectangu- 
lar, sixteen-degree-of -freedom  element  was  used.   Results  show  that, 
for  the  same  order  of  error,  this  element  allows  much  coarser  plate 
divisions  than  a  widely  used  twelve-degree-of-freedom  element.   An 
order  of  magnitude  reduction  in  the  required  order  of  the  assembled 
stiffness  and  stability  coefficient  matrices  results. 
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LIST  OF  SYMBOLS 

a  plate  length  in  x  direction;  polynomial  coefficient 

a  vector  of  displacement  polynomial  coefficients 

A  transformation  matrix 

b  plate  length  in  y  direction 

B  stability  coefficient  matrix 

D  plate  flexural  rigidity,  Eh3/(12(1  -   v  2)) 

E  modulus  of  elasticity 

t  vector  of  external  nodal  forces  and  moments 

G      matrix  of  functions  of  x  and  y,  relating  vectors  a 
and  £ 

h      plate  thickness 

H      matrix  of  functions  of  x  and  y,  relating  vectors  a 
and  r_ 

1  identity  matrix 

k  buckling  coefficient 

K  stiffness  matrix 

K*  modified  stiffness  matrix  from  defining  relation  K  =  DK* 

m  vector  of  functions  of  bending  moments 

M  bending  moment  per  unit  length 

N*  load  matrix 

N      matrix  of  constants  giving  relative  magnitudes  of  types 
of  loads 

p  plate  element  length  in  x  direction 

P  matrix  of  functions  of  Poisson's  ratio 

q  plate  element  width  in  y  direction 

r  vector  of  first  derivatives  of  transverse  displacement 

U  strain  energy  of  bending 

V  total  potential  energy 
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w      transverse  displacement 

W      work  of  in-plane  forces 

x,y    rectangular  coordinates 

z      vector  of  second  derivatives  of  transverse 
displacement 

CL  clamped  edge 

FR  free  edge 

SS  simply  supported  edge 

(5  vector  of  nodal  displacements  and  slopes 

i\  scalar  factor  of  matrix  N* 

'X  eigenvalue 

"0  Poisson's  ratio 

T 

transpose  of  a  matrix 

inverse  of  a  matrix 

in  the  x  direction;  for  moments,  about  a  line 
parallel  to  the  x  axis 

partial  derivative  with  respect  to  x 

pertaining  to  a  plate  element 


1.  INTRODUCTION 
BACKGROUND 
The  Finite  Element  Method  for  Stability  Problems 

The  finite  element  method  provides  a  numerical  means  of  analyzing 
boundary  value  field  problems.  When  applied  to  structural  mechanics, 
the  method  provides  a  means  of  deriving  a  stiffness  matrix  for  an  as- 
sembly of  structural  elements.   For  stability  problems,  a  stability  co- 
efficient matrix  is  also  required.   This  matrix  is  used  in  conjunction 
with  the  stiffness  matrix. 
Development  of  the  Method 

Turner,  Clough,  Martin,  and  Topp  introduced  the  finite  element 

1 
method  in  1956.    In  1963,  Melosh  presented  a  basis  for  derivation  of 

the  stiffness  matrix  for  a  thin  plate  element,  using  three  degrees  of 

2 
freedom  per  node.   Also  in  1963,  Gallagher  and  Padlog  first  used  a 

3 
stability  coefficient  matrix,  for  column  buckling.    In  1966,  Kapur  and 

Hartz  gave  the  results  of  the  use  of  the  stability  coefficient  matrix 

4 
approach  for  plate  buckling  problems.   They  used  the  Melosh  plate  ele- 
ment to  find  the  stiffness  matrix.  While  Kapur  and  Hartz  were  making 

their  study  of  plate  buckling,  in  1965  Bogner,  Fox,  and  Schmit  presented 

5 
a  new  stiffness  matrix  for  a  plate  element.   Four  degrees  of  freedom 

per  node  were  used  for  the  new  element. 

Significance  of  Degrees  of  Freedom 

For  the  Melosh  method  of  deriving  the  stiffness  matrix  of  a  plate 

element,  the  three  degrees  of  freedom  per  node  are  lateral  deflection 

and  slope  in  each  of  two  perpendicular,  in-plane  directions  (w;  w,  ; 

w,  ).   This  assures  compatibility  of  lateral  deflection  of  adjacent  edges 

of  elements.   However,  complete  continuity  of  slope  is  not  achieved.   The 
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method  of  Bogner,  Fox,  and  Schmit  includes  an  additional  degree  of 
freedom  (w,  )  which  does  provide  continuity  of  slope.  For  a  rectangu- 
lar plate  element,  an  element  stiffness  matrix  of  order  sixteen  is  re- 
quired to  represent  four  degrees  of  freedom  at  each  of  four  nodes.   For 
three  degrees  of  freedom  per  node,  the  order  is  twelve.   Hence,  the  six- 
teen-degree-of-freedom  element  stiffness  matrix  is  larger  but  should 
produce  more  accurate  results  because  of  better  continuity  of  slope.* 

THE  PROBLEM  INVESTIGATED 
Objective 

The  objective  of  the  investigation  was  to  determine  the  consequences 
of  the  use  of  the  newer,  sixteen-degree-of-freedom  element  stiffness 
matrix  in  plate  buckling  problems.   Secondarily,  some  previously  unsolved 
cases  of  plate  buckling  were  investigated  to  show  the  versatility  of  the 
finite  element  method. 
Limitations 

The  investigation  was  limited  to  finding  buckling  coefficients  and 
buckling  mode  vectors  for  thin,  homogeneous,  isotropic,  rectangular  plates 
with  free,  simply  supported,  or  clamped  edge  conditions,  and  with  uniform 
compression  or  pure  shear  loads.   Plates  were  subdivided  into  identical 
rectangular  elements. 

JUSTIFICATION  OF  THE  INVESTIGATION 
Advantages  of  the  Finite  Element  Method 


*Other  rectangular  plate  elements  with  twelve  or  sixteen  degrees  of 
freedom  may  be  formed  by  using  different  displacement  functions.   Through- 
out, plate  elements  with  twelve  or  sixteen  degrees  of  freedom  will  mean 
the  elements  discussed  above. 
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Plate  buckling  problems  of  even  moderate  complexity  require  numeri- 
cal procedures .  Of  the  numerical  procedures  applicable,  the  finite  ele- 
ment method  has  some  unique  advantages,  as  follows; 

1.  The  capability  of  analyzing  plates  of  arbitrary  shape,  such 
as  plates  with  holes . 

2.  The  capability  of  analyzing  plates  made  of  materials  with 
orthotropic  properties, 

3.  The  capability  of  introducing  edge  conditions  conveniently, 

especially  the  free  edge. 

4.  The  capability  of  synthesizing  complex  structures  from  plates 
and  other  members,  such  as  beams. 

Advantages  of  Using  Sixt:een-j)eg^ree-qfj^F_reedojn  Element 

Assembled  stiffness  and  stability  coefficient  matrices  may  be  very 
large,  especially  in  structural  problems  which  require  synthesis  of 
several  types  of  members.   Finer  divisions  of  a  plate  require  larger 
assembled  matrices.   Any  procedure  which  would  reduce  the  order  of 
these  matrices,  without  increasing  the  error  of  the  result,  would,  of 
course,  save  computer  time  and  space. 

It  was  anticipated  that  use  of  the  newer  element  would  allow  coarser 
plate  divisions  with  no  increase  of  error.   If  so,  a  significant  reduc- 
tion in  the  size  of  the  assembled  matrices  might  be  possible  in  spite  of 
the  higher  order  of  the  element  matrices.   Hence,  the  investigation  was 
undertaken  to  see  whether  the  newer  element  produces  a  significant  im- 
provement. 

PROCEDURE 
The  equations  for  bending  energy  and  work  of  in-plane  forces  were 
written  in  matrix  form,  making  use  of  the  finite  element  method  with 
the  sixteen-degree-of-freedom  element.   A  computer  program,  with  matrix 
iteration,  was  written  to  generate  and  solve  the  matrix  equation  re- 
sulting from  application  of  conditions  of  minimum  potential  energy. 
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Several  previously  solved  cases  of  plate  buckling  were  run.   Results 
were  compared  with  known  values  and  with  values  obtained  using  the  twelve- 
degree-of -freedom  plate  element.   A  few  previously  unsolved  cases  were 
investigated. 
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2.  FORMULATION 
INTRODUCTION 
In  this  section,  the  formulation  of  the  problem  is  introduced  by 
a  discussion  of  the  energy  method.   The  plate  element,  displacement 
vector,  displacement  polynomials  and  transformation  matrix  are  explained 
as  preliminaries.  With  this  background  establisheds  the  matrix  equa- 
tions for  bending  energy  and  work  of  in-plane  forces  are  derived  for  a 
single  plate  element.   Extension  of  the  method  to  multiple-element 
plates  and  the  method  of  handling  boundary  conditions  are  discussed  next. 
Then,  the  method  of  finding  a  matrix  equation  for  minimum  potential  energy 
is  shown.   Finally,  the  solution  of  the  resulting  characteristic  value 
problem  is  discussed. 

THE  ENERGY  METHOD 

If  a  plate  loaded  by  externals  in-plane  forces  is  bent  a  small 
amount,  the  external  forces  move  because  the  plate  edges  move.   Thus, 
the  forces  do  work.   At  the  same  time,  to  bend  the  plate,  the  strain 
energy  of  bending  must  be  supplied.   If  the  work  of  external  forces 
exceeds  the  bending  energy,  the  plate  buckles.   These  two  energy  terms 
are  equal  in  magnitude  at  the  critical  load.   Or,  the  total  potential 
energy  is  minimum  at  the  critical  load. 

PRELIMINARY  DEFINITIONS 

Plate  Element  and  Displacement  Vector 
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—  X 


The  plate  element  used  is  shown  above.   The  nodal  order  used  is 
indicated  at  the  corners  of  the  element.   The  generalized  displacement 
vector  used  for  a  single  node  is* 


61  =  C  w  w,y  w,K  w,XylT 


(1) 


The  displacement  vector  used  for  a  plate  element  is 


An  underline  indicates  a  matrix.   Capital  letters  are  used  for 
square  or  rectangular  matrices.   Small  letters  are  used  for  column 
vectors.   A  superscript  T  is  used  to  indicate  a  transpose. 
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^o  -  I  tf±(0,0)  j  £(Pj0>  ;  ^0,q)   ;  6A(p,q)3r  (2) 

Displacement  Polynomial  and  Transformation  Matrix 
Transverse  deflection  was  represented  by 

vt-i  1  ayx^y5"*  (3) 

i=i  j=i  J 

By  using  this  polynomial  and  derivatives  of  it,  evaluated  at  the  plate 
element  nodes,  the  generalized  deflection  of  the  nodes  may  be  represented 
by 

do  =  Ag  (4) 

where  A  is  a  16  x  16  transformation  matrix,  and  a  is  a  vector  of  the  poly- 
nomial coefficients.   The  polynomial  coefficients  are  then  given  by 


3  =  AT^o  (5) 


BENDING  ENERGY 
The  integral  expression  for  the  strain  energy  of  bending  is 

U  =   -|JJ  (Mxw,xx  +  Myw/yy  -2.MxywjXy)dxdy 


(6) 


This  expression  may  be  re-stated  in  matrix  form  as 

U  =    -iJjmTz  dxdy  (7) 
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where  m  and  z  are  vectors  defined  by 


nn=[Mx   My  -ZMxy3T 

Z     =    CWJXX     W>yy      W,Xy] 


(8) 


(9) 


Bending  moments  are  related  to  curvatures  and  material  properties  by 

Mx  =  -  D  (wjXx  +  V  W;yy) 
My=~D(W,yy  +  "V  W,xx) 
M*y=  OU-'O)  Wyxy 

Here,  D  is  the  plate  flexural  rigidity,  and  -0  is  Poisson's  ratio. 
Equations  (10)  may  also  be  written  as 


(10) 


Mx 

My 
-2M*y 


=  -D 


i 


0 


W,yy 
W/xy 


(11) 


V    I    0 

0    0  2(1-0) 

Defining  the  matrix  of  functions  of  0   to  be  P  allows  equation  (11) 
to  be  given  in  the  more  compact  form 


m  = -QEz 


(12) 


Combining  equation  (7)  and  equation  (12)  gives  the  new  equation  for 
bending  energy 


U  =   |DjJ(Pz)Tz  dxdy 


(13) 


Since  the  vector  z   is  made  up  of  derivatives  of  the  displacement  poly- 
nomial, it  may  be  related  to  the  polynomial  coefficients  by 

£  =  £5  (14) 

where  G  is  a  3  x  16  matrix  of  functions  of  x  and  y.   Using  equation  (14) 
to  substitute  for  z   in  equation  (13)  gives 


U  =  |  DjJYPGafSg  dx  dy 
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(15) 


Substituting  for  the  vector  a  using  equation  (5) ,  extracting  terms  not 
containing  x  and  y  from  the  integrand,  and  rearranging  gives 


The  portion 


TJ  =   |  d  60T  (A-1)1"^  GT  P  G  dx  dy )  A_1d0 


K„  =  D^VljjVp/Gdxd^A'1 


(16) 


(17) 


is  the  stiffness  matrix  for  the  plate  element.   Finally,  the  bending 
energy  of  the  plate  element,  in  terms  of  the  element  stiffness  matrix, 


is 


U  =  I  &Z  K'o  £< 


(18) 


WORK  OF  IN-PLANE  FORCES 
The  integral  expression  for  the  additional  work  done  by  in-plane 
forces  because  of  lateral  deflection  is 


The  matrix  form  of  this  equation  is 


W-  -|JJrTN*  rdxdy 


(20) 


Here,  r  and  N*  are  defined  by 


r  = 


W,y 


N*= 


Ny  Nxy 
_Nxy  Nx 


(21) 


(22) 
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The  matrix  of  load  parameters  may  be  expressed  as 

N*=r\,N  (23) 

where  n   is  a  scalar  factor,  and  N  is  a  matrix  of  constants  which  give 
the  relative  magnitudes  of  the  types  of  loads.  The  vector  r,  which  con- 
sists of  derivatives  of  the  displacement  polynomial,  may  be  related  to 
the  polynomial  coefficients  by 

r  =  Hd  (24) 

H  is  a  2  x  16  matrix  of  functions  of  x  and  y.   Equations  (20),  (23), 
and  (24)  lead  to 


W  =  -  |»\I[(Ha)TNHgdxdy  (25) 


Continuing  in  the  same  manner  as  was  used  for  bending  energy  formulation 
gives 

W  «  -  £  i\  dj  (A"1)1  (JjHTNM  dx  dy )  A"ld0  (26) 

The  portion 

B0  =  (A"MT  (|J HT  N  H  dx  dy  )  A"1  (27) 

is  the  stability  coefficient  matrix  for  the  plate  element.* 

The  work  of  in-plane  forces,  in  terms  of  the  stability  coefficient  matrix, 


is 


W*    -|  >\6jB0rfo  (28) 


*The  stability  coefficient  matrix  is  sometimes  referred  to  as  the 
geometric  stiffness  matrix. » 
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EXTENSION  TO  MULTIPLE  ELEMENTS 

So  far,  the  stiffness  matrix  and  stability  coefficient  matrix  for 
a  single  plate  element  have  been  derived.   Poor  accuracy  results  when  a 
plate  is  represented  by  a  single  element.   Rather,  the  plate  must  be 
represented  by  an  assembly  of  plate  elements.   Although  the  plate  need 
not  be  divided  into  identical  elements,  identical  elements  were  used 
throughout  this  investigation. 

The  same  method  of  extension  to  multiple  elements  for  the  stiffness 
matrix  applies  to  the  stability  coefficient  matrix.   Only  a  brief  dis- 
cussion will  be  given.   A  more  complete  description  is  contained  in 
references  7  and  8. 

Let  f  be  a  generalized  force  vector  representing  externally  ap- 
plied forces  and  moments  at  nodes  of  a  plate  element.   Then,  f  and  the 
corresponding  nodal  displacement  vector  are  related  by 

fo  =  Ko  do  (29) 

where  K  has  an  order  equal  to  the  number  of  degrees  of  freedom  of  the 
plate  element.  Similarly,  when  several  elements  are  assembled  to  make 
up  a  plate, 

f  =  Ktf  (30) 

where  K  is  the  assembled  stiffness  matrix,  of  higher  order  than  K  be- 
—  — o 

cause  there  are  more  nodes.   Elements  of  K  are  sums  of  elements  taken 
from  the  plate  element  stiffness  matrices. 

BOUNDARY  CONDITIONS 
One  of  the  appealing  characteristics  of  the  finite  element  method  is 
the  simplicity  of  handling  boundary  conditions.   For  example,  no  action 
whatever  is  required  for  a  free  edge.   For  other  edge  conditions,  two  ap- 
proaches are  available.   One  method  is  to  eliminate  appropriate  rows  and 

columns  of  the  stiffness  and  stability  matrices,  reducing  the  order  of  the 
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matrices.   For  ease  of  computer  programming,  a  second  method  was  used. 
Zeros  were  placed  in  appropriate  rows  and  columns  of  the  matrices.   Then 
each  zero  diagonal  element  of  the  stiffness  matrix  was  made  unity,  be- 
cause the  inverse  of  the  matrix  was  needed  later.   Table  I  shows,  for 
edge  nodes,  the  nodal  displacement  components  which  are  zero  for  the 
prescribed  edge  condition.   It  should  be  noted  that  the  constraints  given 
in  Table  I  are  the  only  edge  constraints.   There  is  no  constraint  on  in- 
plane  deflection  for  any  edge  condition;  i.e.,  the  plate  is  "supported" 
or  "clamped",  but  not  "held". 

MINIMUM  POTENTIAL  ENERGY 
Let  K  and  B  now  represent  the  assembled  stiffness  and  stability 
matrices,  altered  by  edge  conditions.   The  total  potential  energy  of 
the  plate  is  the  sum  of  bending  energy  and  work  of  in-plane  forces. 
From  equations  (18)  and  (28), 

V=  U*W  =  46TK5  -  |r\rfTBrf  (31) 

The  critical  condition  for  stability  is  an  equilibrium  configuration 
with  a  set  of  displacements  which  makes  the  total  energy  a  minimum.   This 
requires  that 

Kd  -  r^Brf  *0  (32) 

CHARACTERISTIC  VALUE  PROBLEM 

With  some  alteration  of  form,  equation  (32)  expresses  a  well-known 

9 
characteristic  value  problem. 

(DK*)6  =ry&6  (33) 
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TABLE  I 


BOUNDARY  CONDITIONS 


Edge  Condition 


w 


w». 


w, 


w, 


xy 


Simply  supported  edge 
parallel  to  x  axis 

Simply  supported  edge 
parallel  to  y  axis 

Clamped  edge 

Line  of  symmetry 
parallel  to  x  axis 

Line  of  symmetry 
parallel  to  y  axis 

Line  of  ant i -symmetry 
parallel  to  x  axis 

Line  of  anti-symmetry 
parallel  to  y  axis 


0 
0 
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(D/r{)&6=  B6  (34) 

'X  K*6  =  &6  (35)* 

Here,  ^\   is  an  eigenvalue  parameter  which  may  have  as  many  values  as 
the  order  of  the  matrices.   Since  only  the  lowest  load  is  desired, 
however,  only  the  largest  eigenvalue  needs  to  be  found. 

Most  plate  buckling  coefficients  are  expressed  in  the  form 

k  =  (Ncrii  b*)/(ir*  D)  (36) 

where  b  is  the  plate  width  in  the  y  direction.   If  the  largest  value 
of  X   is  known,  the  buckling  coefficient  may  be  given  in  the  same 
form  by 

k=(rlb2)/(TT2D)=  bVOn-*)  (37) 

SOLUTION  OF  CHARACTERISTIC  VALUE  PROBLEM 
Solving  equation  (35)  for  the  largest  eigenvalue  is  the  only 

problem  remaining.   Matrix  iteration  was  the  method  used,  although 

there  are  other  methods. 

The  first  procedure  tried  was  as  follows: 

1.  Assume  a  starting  trial  vector,  (£  . 

2.  Calculate  an  improved  trial  vector,  6',  using 

6'  =  K-1Bd  (38) 

3.  Use  the  Rayleigh  quotient  to  find  an  estimate  for   "A.  ,  using 

A=  (^TBtf'J/UdyKd')  (39) 

4.  Compute  k,  using 


,-ti 


*This  equation  is  often  given  in  the  form  ((K*)  B->J)d  =  G. 
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k  =  bV(*ir2)  (40) 

5.   Use  the  improved  trial  vector  as  the  next  starting  trial 
vector,  and  iterate  until  successive  values  of  k  agree 
within  a  prescribed  limit. 

The  procedure  did  not  give  good  results  for  shear  loads.   For  pure  shear, 
corresponding  to  each  positive  eigenvalue,  there  is  a  negative  eigen- 
value of  equal  magnitude.   The  Rayleigh  quotient  does  not  give  a  good 
estimate  for   "A.  . 

A  new  procedure  was  used  which  proved  to  be  satisfactory.   A  new 
quotient, 

A2  =((6,)TKd,)/(dTKd)  (41) 

was  used  instead  of  the  Rayleigh  quotient.   Otherwise,  the  procedure 
was  the  same  as  before. 
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3.  RESULTS 
INTRODUCTION 
General 

In  this  section,  results  of  three  kinds  are  given.   First,  compari- 
sons are  made  with  known  solutions  to  show  the  accuracy  of  the  method. 
Second,  comparisons  are  made  with  results  given  by  Kapur  and  Hartz  to 
show  the  effect  of  the  sixteen-degree-of -freedom  element.   Third,  some 
previously  unsolved  cases  of  buckling  are  given  to  show  the  versatility 
of  the  finite  element  method. 
Information  for  Reading  Tables  of  Results 

Wherever  appearing,  the  known  values  of  the  buckling  coefficient 
are  from  the  methods  given  by  Timoshenko  and  Gere,  but  are  not  necessar- 
ily  the  values  stated  by  them.   Some  values  were  computed  to  a  greater 
precision.   In  all  cases,  the  buckling  coefficient  is  given  in  the  form 

k=(NcKt  t?)/(ir*D)  (42) 

ACCURACY  OF  THE  RESULTS 
General 

Table  II  shows  the  effect  of  grid  size  upon  the  error  of  the  buck- 
ling coefficient  for  several  cases.   Table  III  shows  the  effect  of  as- 
pect ratio  upon  error  for  a  single  case.   In  many  cases,  the  error  is 
less  than  one  per  cent  for  rather  coarse  grid  sizes. 
Factors  Affecting  Error 

General  observations  concerning  the  results  are  as  follows: 

1.  Error  is  greater  for  coarser  grids,  as  expected. 

2.  Error  is  greater  for  pure  shear  loading. 

3.  Error  is  greater  for  clamped  edges. 

4.  Error  is  greater  for  narrower  or  longer  plates. 
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TABLE  II 

COMPARISON  OF  BUCKLING  COEFFICIENTS  WITH  KNOWN  VALUES  FOR 
SEVERAL  CASES,  USING  SEVERAL  GRID  SIZES 


Uniaxial  Compression 
Aspect  Ratio  =  1 
Poisson's  Ratio  =  0 
k(knovm)  -  4.00000 


Grid  Size 

k 

% 

Difference 

2x2 

4.01575 

+0.4 

3x3 

4.00324 

+0.08 

4x4 

4.00104 

+0.03 

6x6 

4.00021 

+0.005 

8x8 

4.00007 

+0.002 

Isotropic 

Compression 

Aspect  Ratio 

=  1 

Poisson's 

Ratio  -  0 

k( known)  s 

■  2. 

00000 

Grid  Size 

k 

% 

Difference 

2x2 

2.00791 

+0.4 

4x4 

2.00052 

+0.03 

8x8 

2.00003 

+0.002 

5S 
— x 


I  ksl  I 


yfsso" 


ss 


CL 


it  fcki  i 


T~FTT 


CL 


Isotropic  Compression 
Aspect  Ratio  =  1 
Poisson's  Ratio  =  0 
5.30<  k(known)  <C  5.33--lower 
bound  used  for  comparison. 


Grid  Size 

k 

% 

Difference 

2x2 

5.36223 

+1.2 

3x3 

5.36599 

+1.2 

4x4 

5.32713 

+0.5 

6x6 

5.30900 

+0.2 

8x8 

5.30544 

+0.1 
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TABLE  II  (continued) 


SS 


— -ss — - 

1 

SS 

Pure  Shear 
Aspect  Ratio 
Poisson's  Rat 
k (known)  =  9. 

Grid  Size 

2x2 
3x3 

4x4 

=  1 

:io  =  0 
34 

k 

%   Difference 

10.01581 
9.57700 
9.41822 

+7.2 

— ss- — 

+2.5 
+0.8 

SS 


CL 


SS 


FR 


Uniaxial 

Compression 

Aspect  Ratio 

=  1 

Poisson's 

Ratio  =  0.25 

k (known) 

=  1. 

70 

Grid  Size 

k 

%   Difference 

2x2 

1.71399 

+0.8 

4x4 

1.69942 

* 

8x4 

1.69895 

* 

*Zero  for  precision  of  known  value. 


— ^CL — »- 

1 

CL 

f 

Pure  Shear 
Aspect  Rat 
Poisson's 
k(known)  ■ 

Grid  Size 

2x2 
3x3 

4x4 

io  =  1 
Ratio  =  0 
14.71 

k 

CL 

%  Difference 

1 

23.26396 
16.04630 
15.04271 

+58.2 

—  CL" 

+9.1 
+2.3 
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TABLE  II  (Continued) 


ss 

SS 

Uniaxial  ( 
Aspect  Rat 
Poisson's 
k  (known ) 

Grid  Size 

2x2 

4x4 
8x4 

Compression 
:io  =  1 
Ratio  »  0.25 
=  1.43418 

k 

ss 

7.  Difference 

1.44324 
1.43483 
1.43435 

+0.6 

+0.05 

+0.01 

FR 

—  CL-^ 

f 

SS 

[ 

Pure  Sheai 
Aspect  Rat 
Poisson's 
k (known)  ■ 

Grid  Size 

2x2 

4x4 

:io  =  1 
Ratio  =  0 
■  12.28 

k 

1 

ss 

1 

7»  Difference 

17.13879 
12.82566 

+39.6 

"-   LL  " 

+4.4 

SS 


CL 


Uniaxial  Compression 
Aspect  Ratio  =  1 
Poisson's  Ratio  =  0 
k (known)  =6.74 


CL 


SS 


Grid  Size 

k 

7o 

Difference 

2x2 

6.82334 

+1.2 

4x4 

6.78242 

+0.6 

8x8 

6.74583 

+0.09 

Uniaxial  Comp 

ression 

Aspect  Rat 

:io 

=  1 

Poisson's 

Rat 

io  =  0 

k (known)  = 

=  7. 

69 

Grid  Size 

k 

7» 

Difference 

2x2 

8.79292 

+14.3 

4x4 

7.74665 

+0.7 

8x8 

7.69541 

+0.07 

CL 


SS 


ss 


CL 
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TABLE  III 


COMPARISON  OF  BUCKLING  COEFFICIENTS  WITH  KNOWN  VALUES 
FOR  ONE  CASE,  USING  SEVERAL  ASPECT  RATIOS 


*  -r 


ss 

X 

Uniaxial 
Grid  Size 

Comp 
,  4 

ression 
x  4 

_,.!    c\            **_ 

-*i 

Poisson's 

Rat 

io  ■  0 

ss 

b 

SS 

Aspect  ratios 

for  buckling  mode 

-*■ 

* 

*— 

transition: 

2*,  3%, 

4*, 

etc. 

ss 

T 

a/b 
0.2 

k 

k (known) 
27.04000 

%  Error 
+0.048 

No.  Half  Waves 

27.05296 

0.4 

8.41330 

8.41000 

+0.039 

0.6 

5.13940 

5.13778 

+0.032 

0.8 

4.20364 

4.20250 

+0.027 

1.0 

4.00104 

4.00000 

+0.026 

1.2 

4.13556 

4.13434 

+0.030 

1.4 

4.47150 

4.47022 

+0.029 

2.0 

4.00770 

4.00000 

+0.193 

2 

** 

3.0 

4.02924 

4.00000 

+0.731 

3 
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When  the  buckled  shapes  of  the  plates  are  considered,  all  of  these 
observations  may  be  further  generalized  to  a  single  one.   The  accuracy 
of  representing  the  buckled  surface  is  less  for  greater  curvature  gradi- 
ents.  To  keep  error  low,  finer  grids  must  be  used  for  more  complex 
buckled  shapes. 
Sources  of  Error 

Discretization  Error  Most  numerical  methods  replace  continuous 
equations  with  discrete  approximations.   The  finite  element  method, 
however,  uses  continuous  equations  but  a  discrete  approximation  of  the 
continuous  structure  being  analyzed.   The  finite  element  method  has  in- 
herent discretization  error.   Since  the  finite  element  approximation 
of  a  real  plate  is  stiffer  than  the  real  plate,  the  discretization 
error  of  the  buckling  coefficient  is  never  negative.   In  addition,  the 
discretization  error  decreases  montonically  as  the  plate  is  further 

subdivided,  provided  that  previously  existing  nodes  are  also  nodes  in 

2  5 
the  finer  grid.  '   For  errors  shown  in  Tables  II,  III,  and  IV,  most 

error  is  believed  to  be  discretization  error. 

Manipulation  Error  Matrix  iteration  produced  monotonically  con- 
verging values  of  the  buckling  coefficient.   Further,  as  a  consequence 
of  the  matrix  iteration  method  used,  convergence  was  always  from  an 
initially  higher  value.   Hence,  matrix  iteration  error  of  the  buckling 
coefficient  is  also  positive.   A  convergence  criterion  of  .001  per  cent 
difference  between  successive  values  of  the  buckling  coefficient  was 
used.   Consequently,  the  error  of  the  results  cannot  be  much  less  than 
.001  per  cent.   Round-off  errors  are  believed  to  be  insignificant  in 
comparison  with  other  errors,  since  double  precision  programming  was 
used. 
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Negative  Error 

Because  the  two  principal  sources  of  error--discretization  and 
convergence --both  produce  positive  error  of  the  buckling  coefficient, 
the  error  cannot  be  negative.   Although  some  errors  initially  appeared 
to  be  negative,  in  each  case  encountered  the  cause  was  found  to  be  lack 
of  precision  of  the  "known"  value  of  the  buckling  coefficient. 
Accuracy  of  Buckling  Mode  Vectors 

Although  the  buckled  shape  of  the  plate  is  of  secondary  interest, 
the  buckling  mode  vector  was  available  as  a  consequence  of  the  method  used. 
The  vector  provided  another  check  on  the  accuracy  of  the  method.   For 
example,  the  nodal  displacements  were  compared  with  values  computed  from 
the  known  buckled  shape  for  the  uniaxially  compressed,  simply  supported 
plate  of  aspect  ratio  one  for  a  4  x  4  grid.   The  displacements  were  in- 
ternally consistent  to  a  minimum  of  three  significant  figures. 

EFFECT  OF  SIXTEEN -DEGREE -OF -FREEDOM  ELEMENT 

Table  IV  shows  comparisons  of  results  with  those  given  by  Kapur  and 

4 
Hartz  for  two  cases.   An  improvement  in  accuracy  is  apparent.  More 

important,  for  the  same  order  of  magnitude  of  the  error,  fewer  elements 
were  required. 

To  show  what  this  improvement  means,  consider  the  case  of  the  simply 
supported  plate  compressed  in  one  direction.   Kapur  and  Hartz  report  an 
error  of  -0.58  per  cent  when  using  a  12  x  12  grid.   Assuming  that  no 
symmetry  and  no  compaction  of  matrices  is  used,  and  not  considering  re- 
duction of  the  matrix  order  when  boundary  conditions  are  applied,  their 
formulation  requires  a  507  x  507  assembled  stiffness  matrix  to  represent 
three  degrees  of  freedom  at  each  of  169  nodes.   A  stability  coefficient 
matrix  of  the  same  order  is  also  required,  of  course.   Using  one  more 
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TABLE  IV 

COMPARISON  OF  BUCKLING  COEFFICIENTS  WITH  VALUES  OBTAINED 
USING  TWELVE -DEGREE -OF -FREEDOM  ELEMENT3 


2x2 
3x3 
4x4 
6x6 
8x8 
10  x  10 
12  x  12 


Uniaxial  Compression 
Aspect  Ratio  =  1 
Poisson's  Ratio  =  0 
k (known)  =  4.00000 


12  Degrees  of  Freedom 


k 

7.  Error 

***** 

**** 

3.645 

-8.88 

3.77 

-5.75 

3.887 

-2.80 

3.933 

-1.68 

3.96 

-1.00 

3.977 

-0.58 

16  Degrees  of  Freedom 


k 

7»  Error 

4.01575 

+0.39 

4.00324 

+0.08 

4.00104 

+0.03 

4.00021 

+0.005 

4.00007 

+0.002 

******* 

'fc'fc'k'k'k'k 

******* 

jcic'k'k'k'k 

CL 


LioJJ 


*  tc4  T 

Grid  Size 


CL 


2 
3 
4 
6 
8 
10  x 


2 
3 
4 
6 
8 
10 


Isotropic  Compression 
Aspect  Ratio  =   1 
Poisson's  Ratio  =  0       . 
5.30<  k(known)<  5.33 


12  Degrees   of  Freedom 


k 

7»  Diff. 

4.901 

-7.53 

***** 

***** 

4.975 

-6.13 

5.078 

-4.19 

5.160 

-2.64 

5.217 

-1.57 

16  Degrees   of  Freedom 


7o  Diff, 


5.36223 
5.36599 
5.32713 
5.30900 
5.30544 
******* 


+1.17 
+1.24 
+0.51 
+0.17 
+0.10 

-t . . #  .»  . ■ . .» 


Data  from  Kapur  and  Hartz,  reference  4. 

The  lower  bound  was  used  for  computing  per  cent  differences. 
Since  Kapur  and  Hartz  had  used  an  average  of  the  upper  and  lower  bounds, 
the  per  cent  differences  given  by  them  were  re-computed. 
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degree  of  freedom  per  node,  the  error  is  0.39  per  cent  when  using  a 
2x2  grid.   For  the  same  assumptions,  the  stiffness  matrix  size  is  36  x 
36  for  four  degrees  of  freedom  at  each  of  nine  nodes.   In  this  example, 
the  superior  accuracy  of  the  newer  plate  element  and  the  consequent 
savings  in  computer  time  and  space  are  unequivocally  apparent. 

PREVIOUSLY  UNSOLVED  PROBLEMS 

Table  V  shows  the  results  of  investigation  of  two  previously  un- 
solved plate  buckling  problems. 

An  edge  which  has  more  than  one  edge  condition  along  its  length 
is  difficult  to  analyze  by  other  means.   Although  the  computer  program 
used  for  other  cases  had  to  be  altered  very  slightly  to  handle  such  an 
edge,  there  was  no  particular  difficulty  in  using  the  finite  element 
method. 
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TABLE  V 
PLATE  BUCKLING  PROBLEMS  NOT  PREVIOUSLY  INVESTIGATED 


FR 

FR 

X 

X 

X 

1 

Ti 

i 

i 

14 

5' 

FR 

X 

il 

1G> 
Zi 

Tit 

1« 

il5 

l 

i 

1 

125 

19 
1l4 

!Z4 

-a 

15 
25 

•am ■ 

FR 

X 

FR< 

A 

FR       FR        X        X        X 

Uniaxial  compression,  aspect  ratio  =  2,  Poisson's  ratio  =0.3 
Nodes  marked  X  are  either  all  simply  supported  or  all  clamped. 


Node 


Node 


w 


HALF  FREE,  HALF  SIMPLY  SUPPORTED 
k  =  0.179 
Relative  Transverse  Displacements  of  Nodes 
Node     w     Node      w     Node 


w 


w 


HALF  FREE,  HALF  CLAMPED 
k  =  0.208 
Relative  Transverse  Displacements  of  Nodes 
Node     w     Node       w    Node     w 


Node   w 


1 

-.3562 

2 

-.1207 

3 

0 

4 

0 

5 

0 

6 

-.3662 

7 

-.1354 

8 

-.0139 

9 

+.0019 

10 

0 

11 

-.3699 

12 

-.1408 

13 

-.0182 

14 

+.0027 

15 

0 

16 

-.3662 

17 

-.1354 

18 

-.0139 

19 

+.0019 

20 

0 

21 

-.3562 

22 

-.1207 

23 

0 

24 

0 

25 

0 

Node   w 


1 

-.3174 

2 

-.0956 

3 

0 

4 

0 

5 

0 

6 

-.3257 

7 

-.1073 

8 

-.0034 

9 

+.0013 

10 

0 

11 

-.3287 

12 

-.1111 

13 

-.0053 

14 

+.0026 

15 

0 

16 

-.3257 

17 

-.1073 

18 

-.0034 

19 

+.0013 

20 

0 

21 

-.3174 

22 

-.0956 

23 

0 

24 

0 

25 

0 
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4.  THE  COMPUTER  PROGRAM 
GENERAL 
A  sample  computer  run  is  given  in  Appendix  I.   The  program  listing 
is  largely  self-explanatory.   The  program  is  in  FORTRAN  IV  with  double 
precision.   The  IBM  360  computer  at  the  Naval  Postgraduate  School  was 
used.   This  computer  has  a  core  memory  capacity  of  512,000  bytes.   The 
HASP  compiler  used  required  about  55,000  bytes. 

CAUTIONS 
Program  Efficiency 

The  program  was  written  as  an  aid  in  producing  this  thesis  only, 
by  a  relatively  inexperienced  programmer.   Paradoxically,  one  purpose 
of  the  thesis  is  to  show  how  computer  time  and  space  can  be  saved,  yet 
both  were  often  sacrificed  for  ease  of  programming. 
Convergence  in  Matrix  Iteration 

Initially,  a  trial  vector  with  a  high  content  of  some  mode  shapes 
was  used  to  start  the  matrix  iteration  process.   This  sometimes  caused 
convergence  to  an  incorrect  mode.   So  that  the  trial  vector  would  not 
be  prejudiced  toward  particular  modes,  a  new  trial  vector  consisting  of 
elements  taken  from  a  table  of  random  numbers  was  used.   In  all  cases 
run  using  the  new  trial  vector,  where  the  proper  buckling  mode  was 
known,  convergence  was  to  the  correct  mode. 

Convergence  was  very  slow  for  cases  with  the  two  largest  eigen- 
values nearly  equal.   Although  methods  of  accelerating  convergence  were 
tried,  no  generally  acceptable  method  was  found.   The  program  is  not 
recommended  for  such  cases. 
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5.  SUMMARY 
CONCLUSIONS 
For  plate  buckling  problems  using  the  finite  element  method, 
the  sixteen-degree-of-freedom  element  is  far  more  accurate  than  the 
twelve-degree-of-freedom  element.   Not  only  is  the  accuracy  superior, 
but  it  is  sufficiently  greater  so  that  sizes  of  assembled  stiffness 
and  stability  coefficient  matrices  may  be  greatly  reduced.   Far  fewer 
plate  elements  are  needed  for  equivalent  accuracy.   Computer  time  and 
space  requirements  are  significantly  reduced. 

PROBLEMS  NEEDING  FURTHER  STUDY 
The  following  problems  were  beyond  the  scope  of  this  investigation: 

1.  An  attempt  should  be  made  to  find  a  general  rule  for  sub- 

dividing the  plate  for  minimum  error.   Elements  of  un- 
equal sizes,  or  more  divisions  in  one  direction  might  be 
used,  for  example. 

2.  A  systematic  investigation  of  many  of  the  remaining  unsolved 

cases  of  plate  buckling  should  be  made. 

3.  Further  study  of  methods  of  determining  the  eigenvalue 

should  be  made  in  order  to  find  a  more  rapid  method. 

4.  A  more  general  and  more  efficient  computer  program 

should  be  devised  for  industrial  use. 
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APPENDIX  I 
SAMPLE  COMPUTER  RUN 
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I       Ck!5HfiES?RHfAPL^El5niWoBnBi?Hc5LigGx056RlfI^VLt?i^ft^?* 

C      MATRIX. 

C 

C      MAXIMUM  PLATE  SIZE  4  ELEMENTS  X  4  ELEMENTS  OR  EQUIVALENT. 

C 

1  IMPLICIT  REAL  *  8  (A-H.O-Z) 

2  DIMENSION  POLYK  16),P0LY2(  16)  , POLY3 ( 16  I , POL Y4<  16 ) , POL Y5 ( 16 ) , 
1PPLY6I 16),A(16,16),AINV(16,16),SELMT(16,16),P(3,3),GCOEF(3,16), 
2GXEXP(3,16),GYEXP(3,16),GINT(16,16) ,S(100,100) ,SINV( 100,100), 
3SELMT (16,16)  .HCtlEF  (  2, 16)  ,  HXEXP  (  2  ,  16  )  .  HYEXP<  2  ,  16  )  ,  H  I  NT  (  16,  16  )  , 
4B<100,100) ,C( 2,2), W( 100) ,WNEXT< 100) .STIFF (100, 100), WS AVE (100) 

201  EQUIVALENCE  ( A( 1 ) ,G I  NT ( 1 ) ,Hl NT ( 1 ) , R( 1 ) , S ( 1 ) ) , ( GCOFF ( 1 ) , HCOEF ( 1 ) , 
1SELMT( 1>,BELMT(1),WSAVE( 1) >,(GXEXP( I) ,HXEXP(1)  ) ,(GYEXP(  1), 
2HYFXP( 1) ) ,(W(1), POLYK 1) ),(WNEXT(1) , POL Y2 ( 1 ) ) , ( P ( 1 ) , ST  I FF ( 1 ) ) 

C 

C  INPUT  FROM  FIRST  DATA  CARD.  ASPFCT  RATIO  (AR),  POISSON'S  RATIO 

C  (POISS),  NO.  ELEMENTS  IN  X  DIRECTION  (MELMT),  NO.  ELEMENTS  IN  Y 

C  DIRECTION  (NELMT),  BOUNDARY  CONDITIONS  FOR  Y=0,  Y=B,  X  =  0,  X=A 

C  ( IBND,JBND,KBND,LBND) ,  RELATIVE  SIZES  OF  LOADS  NX,  NY,  NXY  (CI, 

C  C2,C3).  FORMAT  2010.0,613,3010.0. 

C 

C  FOR  BOUNDARY  CONDITIONS,  1=FREE,  2=SIMPLY  SUPPORTED,  3=CLAMPE0, 

C  ASYMMETRY,  5=  ANTI -S  YMME  TRY.  4  AND  5  MAY  BF  USED  ONLY  FflR  JBND 

C  AND  LBND.  NX  AND  NY  MUST  NOT  BE  NEGATIVE  (COMPRESSION  ONLY).  NX 

C  AND  NY  MUST  BE  ZERO  IF  NXY  IS  OTHER  THAN  ZERO  (PURE  SHFAR  ONLY). 

C 

196  READ  (5.197)  AR , PO I SS, MEL MT , NELMT , I BNO, JRND,KBND , LBND, CI , C2 , C3 

197  FORMAT  (2010.0  ,613,3010.0) 
C 

C      WRITE  INPUT  DATA 
C 

198  WRITE  (6,199)  AR , PHI SS , MELMT , NFL MT 

199  FORMAT  ( 1H 1 ,  1 3HA SP ECT  RATI 0=, F8. 5, 3X , 14HP0I SSON  RAT  I  0= , F4 . 2, 3X , 
URHNO.  ELEMENTS  IN  X= , 1 1 , 3 X, 18HN0.  ELEMENTS  IN  Y=,I1) 

202  WRITE  (6,203) 

203  FORMAT  ( 1H0, 69HB0UNDAR Y  CONDITIONS,  1=FREE,  2=PINNED,  3=CLAMPF0,  4 
1=SYM.,  5=ANTI-SYM.) 

WRITE  (6,204)  I  BND , JBND,KBND , L RND 

204  FORMAT  ( 1H0, 18HY=0 , EDGE  CONDI T ION, I  2, 3X , 1 8HY  =  B , FDGE  CONDI T ION, I  2, 
13X,18HX=0,EDGE  C ONDl T ION , I  2 , 3X , 1 8HX  =  A, EDGE  CONDI T I  ON, I  2 ) 

WRITE    (6,3)    C1.C2.C3 

3  FORMAT     ( 1H0.24HRELATIVE    LOAD     SIZES       NX= , F5 . 2 , 2X , 3HNY=, F5. ? , 2X , 
14HNXY=,F5.2) 

C 

C      THIS  PORTION  OF  THE  PROGRAM  PRODUCES  THE  RECTANGULAR  PLATE  ELEMFNT 

C      STIFFNESS  MATRIX. 

C 

205  EM=MELMT 

206  EN=NELMT 
EX=1./EM 
EY=l./( AR*EN) 

IF  ( (LBN0.EQ.4) .OR. (LBN0.EQ.5)  )  EX=  1 .  / ( 2 . *FM) 

IF  ( ( JBND.EQ.4) .OR. ( JBND.E0.5) )  EY=1 . / ( 2 . *EN* AR  ) 

C      FSTABLISH  ARRAYS  FOR  EXPONENTS  AND  MULTIPLIERS  FOR  TFRMS  OF 
C      POLYNOMIAL  W  AND  ITS  DERIVATIVES 

c 

4  K=l 

5  DO  13  1=1,4 

6  DO  13  J=l,4 

7  P0LY1(K)=I-1 

8  P0LY2(K)=J-1 

9  POLY3(K)=I-2 

10  P0LY4(K)=J-2 

11  P0LY5(K)=I-3 

12  P0LY6(K)=J-3 

IF  (POLYl(K).LT.O. )  P0LY1(K)=0. 

IF  <P0LY2(K) .LT.O. )  P0LY2(K|=0. 

IF  (PPLY3(K)  .LT.O. >  P0LY3(K)=0. 

IF  ( POL Y4(K)  .LT.O. )  PnLY4(K)=0. 

IF  (P0LY5(K)  .LT.O. >  POLY5(K)=0. 

IF  (POL Y6(K) .LT.O. )  P0LY6(K)=0. 

13  K  =  K«-1 

C      ESTABLISH  TRANSFORMATION  MATRIX  A. 


C 


14  DO  23  1=1,13,4 

15  DO  28  J  =  l,  16 

16  X=0. 

17  IF  ( ( I.GT.4.AND. I ,LT.9).nR.( I.GT.12) )  X=FX 

18  Y=0. 

19  IF  (I.GT.8)  Y=EY 

20  IF  (  (X.EQ.O.  LAND. (POLYK  J). EQ.O.)  )  GO  TO  23 

21  TERM1=X**P0LY1( J) 

22  GO  TO  24 

23  TERM1=1. 

24  IF  (  ( Y.EQ.O.  )  .AND. (P0LY2( J) .EO.O.)  )  GO  TO  27 

25  TERM2=Y**POLY2( J) 

26  GO  TO  2R 

27  TERM2=i. 

28  A( l,J)=TERMl*TERM2 

29  DO  43  1=2, 14,4 

30  DO  43  J=l,16 

31  X=0. 

32  IF  ( ( I. GT. 4. AND. I .LT.9).nR.( I.GT.12) )  X=FX 

33  Y=0. 

34  IF  (  I  .GT.R  )  Y=EY 

35  IF  ( (X. EO.O. ). AND. (POLYK J). EO.O.) )  GO  TO  38 

36  TERM1  =  X**P0LYK  J) 

37  GO  TO  39 

38  TFRM1=1. 

39  IF  ( (Y.EQ.O. ). AND. (P0LY4( J). EO.O.) )  GO  TO  42 

40  TERM2=Y**P0LY4( J) 

41  GO  TO  43 

42  TERM2=1. 

43  A(I,J)=P0LY?(J)*TERM1*TFRM2*FY 
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c 

c 

c 

36 

87 

88 

89 

90 

91 

Si 

94 

95 

96 

97 

98 

99 

100 

101 

102 

103 

104 

105 

106 

107 

108 

109 

110 

111 

112 

113 

44  DO  58  1=3,15,4 

45  00  58  J=l,16 

46  X=0. 

47  IF  ( ( l.GT.4.AMD.I.LT.9).nR.( I.GT.12) )  X=FX 

48  Y=0. 

49  IF  ( l.GT.8)  Y=EY 

50  IF  (  I  X.EQ.O.  LAND.  (  POL  Y3  (  J  |  .  EO.O.  )  )  00  TO  53 

51  TERM1=X**P0LY3I J) 

52  GO  TO  54 

53  TERM1=1. 

54  IF  ( ( Y.EO.O. ).ANO. (P0LY2I J) .FO.O.)  )  GO  TO  57 

55  TFRM2=Y**POLY2I.I) 

56  GO  TO  58 

57  TERM2=1. 

58  A(I,J»=P0LY1(J)*TERMI*TFRM2*FX 

59  DO  73  1=4,16,4 

60  00  73  J=l,16 

61  X=0. 

62  IF  ( I  I. GT. 4. AND. I. IT. 9). OR. I  I.GT.12))  X  =  FX 

63  Y=0. 

64  IF  ( I.GT.8)  Y=FY 

65  IF  ( (X.EQ.O.  )  .AND. ( POL Y3( J I.EQ.O. I  )  GO  TO  68 

66  TERM1=X**P0LY3( J) 

67  GO  TO  69 

68  TFRMl=l . 

69  IF  ( IY. FO.O. LAND. (POL Y4( J). EO.O. ) )  GO  TO  72 

70  TERM2=Y**POLY4( j) 

71  GO  TO  73 

72  TERM2=i. 

73  A  I  I, J)=P0LY1( JI*POI  Y?( J ) *TERM1 *TERM2*EX*EY 

FIND  INVERSF  OF  MATRIX  A  (AINV).  GAUSS-JORDAN  METHOD. 

DO  90  1=1,16 

DO  90  J  =  l,  16 

AINVI I,J)=0. 

IF  ( I  .EO.J)  AINVI I , J)  =  l . 

CONTINUE 

EPS=1.D-10 

DO  113  1  =  1  ,  1  A. 

K=I 

IF  (1-16)  95,102,95 

IF  (A(I.I)-EPSI  96,97,102 

IF  (-A(  1,1  l-FPS)  97,97,102 

K=K+1 

DO  100  J=l,16 

AINV( I, J)=AINV(I,J)+AINV(K,J) 

A(  I,J)=A(  I  ,J)+A(K, J) 

GO  TO  95 

DIV=A( I ,1 ) 

DO  105  J=1,16 

AINVI I,J)=AINV< I ,J)/DI V 

A(  I.J)  =  A( I  ,J)/DI V 

DO  113  L=l,16 

OELT  =  ML.I  ) 

IF  (DABS(DELT)-EPS)  113,113,109 

IF  <L-l  )  110,113,110 

DO  112  J=l,16 

AINV(L,J)=AINV(L.J)-AINVII ,J)*OELT 

A(L,J)=A(L,J)-A(I,J)*DELT 

CONTINUE 

00  843  1=1,16 

DO  843  J=l,16 

IF  (OABSI  1INVI  I,  J)  ).LT. 1.0-10)  AINV(I,J)=0. 

CONTINUE 

ESTABLISH  MATRIX  P,  THE  "ATRIX  OF  FUNCTIONS  OF  POISSON'S  RATIO 

PU,1)  =  1. 

P(  1,2)=P0ISS 

P(  1,3)  =  0. 

P(2, 1)=P0ISS 

P(2,2)=l. 

P( 2,3)=0. 

P(3,l)=0. 

P( 3,2)=0. 

P( 3,3)=2.*<  l.-POISS) 


843 


126 
127 
128 
129 
130 
131 
132 
133 


ESTABLISH  MATRIX  G 
134  DO  143  1=1,16 


135 
136 
137 
138 
139 
140 
141 
142 
143 


GCOEFd.I  ) 
GC0EF(2,I ) 
GCOEFI 3, I ) 
GXEXPI 1,1 ) 
GXEXPI2.I ) 
GXEXPI 3, I  ) 
GYEXPI 1,1 ) 
GYFXP(2,I) 
GYEXPI 3,1 ) 


PPLYK  I 

P0LY2II 

POLYK  I 

PPI.Y5I  I 

PCLYK  I 

=PCLY3( I 

=PPLY2( I 

=P0LY6( I 

=P0LY4( I 


*PPLY3<  I  ) 
*P0LY4< I  ) 
♦P0LY2I  I  ) 


THE  FOLLOWING  PART  OF  THE  PROHRAM  COMPUTFS  INTEGRA!  OVFP  X  FRoy 

0  TO  EX,  INTEGRAL  OVER  Y  FROM  0  TP  FY,  OF  T,  TRANSPOSE  *  P  *  G.  THF 

RESULT  IS  NAMED  GINT. 

144  DP  162  1=1,16 

145  DO  162  J=l,16 

146  SUM1=0. 

147  DO  151  K=l,3 

148  SUM2=0. 

149  DP  150  L=1.3 

150  SUM2  =  SUM2*GrPEF(K,  I  )  *P  (  K  ,  I  )  *GCOEF  (  L  ,  J  )  *  (  1  .  /  I  (  GXFX"  (  K  ,  I  )  +OXFXP  ( I.  ,  .1 1 
1*1  .  )*(GYFXP(K,  I)*GYEXP(L  ,  J)*l.  )  )  )*(  FY**(GYPXP(K,  I  )  *GYPXP(L,  J  )«-l.  )  ) 
2*(FX**(GXEXP(K,I )*GXFXP(L ,J)+1. ) ) 

151  SUM1=SUM1+SUM2 
162  GINT( I, J)=SUMl 
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C      FIND  ELEMENT  STIFFNESS  MATRIX  (SELMT).   SELMT=AINV  TR ANSPOSF*G INT* 

C      A I NV 

C 

175  00  183  1=1,16 

176  00  183  J=l,16 

\U  88Mf83-K-i.i6 

179  SUM2=0. 

180  00  181  L=l,16 

181  SUM2=SUM2*GINT(K,L)*AINV(L,J) 

182  SUM1=SUM1*AINV(K,I»*SUM2 

183  SELMT(I,J)=SUM1 
C 

C      GENERATE  ASSEMBLED  STIFFNESS  MATRIX  S.  THE  FIRST  ELEMENTS  OF  A  X  A 
C       ARRAYS  ARE  ASSIGNED  THEIR  PROPER  LOCATIONS.  THE  RFMAINDER  OF  EACH 
C      ARRAY  IS  GENERATED  FROM  THE  FIRST  ELEMENT. 
C 

200  INDEX=(MELMT*1)*(NELMT*1)*4 
11=1 

209  DP  211  1  =  1, INDEX 

210  00  211  J=l, INDEX 

211  S< I,J)=0. 

212  M=MELMT+1 
N=NELMT+1 

214  11=0 

215  DO  372  1  =  1,  M 

216  DO  372  J=1,N 

217  JJ=0 

218  DO  371  K=1,M 

219  00  371  L=1,N 

2  20  IF  ((K.GT.(I*l)).OR.(K.LT.(I-l)).OR.(L.GT.(J+l)).OR.(L.LT.(J-l))) 
1G0  TO  371 

221  IF  (K-I)  222,227,228 

222  IF  (L-J)  233,238.245 

233  IF  (  (  I.EU.D.OR.U.EO.  1)  )  GO  TO  371 

GO  TO  224 

238  IF  (  I.EO.l  )  GO  TO  371 

IF  ( J.EO.l)  GO  TO  230 

IF  (J.EQ.N)  GO  TO  235 

GO  TO  240 

245  IF  ( (  I.EO.l)  .0R.( J.EO.N) >  GO  TO  371 
GO  TO  248 

227  IF  (L-J)  246,251,252 

246  IF  (J.EO.l)  GO  TO  371 
IF  (  I.EO.l)  GO  TO  254 
IF  ( I.EO.M)  GO  TO  259 
GO  TO  264 

251  IF  (  I.EO.l)  GO  TO  257 
IF  (  I.EO.M)  GO  TO  262 
IF  (J.EO.N)  GO  TO  298 
IF  ( J.EO.l)  GO  TO  305 
GO  TO  319 

262  IF  (J.EO.N)  GO  TO  288 

IF  (J.EO.l)  GO  TO  293 

GO  TO  312 

257  IF  (J.EO.l)  GO  TH  271 

IF  (J.EO.N)  GO  TO  27ft 

GO  TO  281 

252  IF  (J.EO.N)  GO  TO  371 
IF  (  I.EO.l)  GO  TO  394 
IF  (  I.EO.M)  GO  TO  397 
GO  TO  400 

228  IF  (L-J)  269,274,279 

269  IF  (( I.EO.M) .OR. (J.EO.l) )  GO  TO  371 

GO  TO  405 
274  IF  ( I  .EO.M)  GO  TO  371 

IF  (J.EO.l)  GO  TO  408 

IF  (J.EO.N)  GO  TO  411 

GO  TO  414 
279  IF  (( I.EO.M) .OR. (J.EQ.N) )  GO  TO  371 

GO  TO  419 

224  111=13 

225  JJJ=1 

226  GO  TO  364 

230  111=5 

231  JJJ=1 

2  32  GO  TO  364 
235  111=13 

237  GO  TO  364 

240  111=5 

241  JJJ=l 

242  KKK=13 

243  LLL=9 

244  GO  TO  354 

248  111=5 

249  JJJ=9 

250  GO  TO  364 

254  111=9 

255  JJJ=1 

256  GO  TO  364 

259  111=13 

260  JJJ=5 

261  GO  TO  364 

264  111=9 

265  JJJ=1 

266  KKK=13 

267  LLL=5 

268  GO  TO  354 

271  111=1 

272  JJJ=1 

273  GO  TO  364 

276  111=9 

277  JJJ=9 

278  GO  TO  364 

281  111=1 

282  JJJ=1 

283  KKK=9 
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284 

LLL=9 

285 

GO  TO  354 

288 

111=13 

289 

JJ.I=13 

290 

GO  TO  364 

293 

m=5 

294 

JJJ=5 

295 

GO  TO  364 

298 

m=9 

299 

JJJ=9 

300 

KKK=13 

301 

LLL=13 

302 

GO  TO  354 

305 

111=1 

306 

JJJ=1 

307 

KKK=5 

308 

LLL  =  5 

309 

GO  TO  3  54 

312 

111  =  5 

313 

JJJ=5 

314 

KKK=13 

315 

LLL=13 

316 

GO  TO  354 

319 

!II  =  1 

320 

JJJ=1 

321 

KKK=5 

322 

LLL=5 

323 

MMM=9 

324 

NNN=9 

325 

1111=13 

326 

JJJJ=13 

327 

GO  TO  340 

394 

m  =  i 

395 

JJJ  =  9 

396 

GO  TO  364 

397 

111  =  5 

398 

JJJ=13 

399 

GO  TO  364 

400 

111=1 

401 

JJJ=9 

402 

KKK=5 

403 

LLL=13 

404 

GO  TO  354 

m  MJ:3 

407 

GO  TO  364 

408 

II  1=1 

409 

JJJ  =  5 

410 

GO  TO  364 

411 

111  =  9 

412 

JJJ=13 

413 

GO  TO  364 

414 

111=1 

415 

JJJ=5 

416 

KKK  =  9 

417 

LLL=13 

418 

GO  TO  354 

419 

111  =  1 

420 

JJJ=13 

421 

GO  TO  364 

340 

00  352  KK=1,4 

341 

MM=I I*KK 

342 

DO  348  LL=1,4 

343 

NN=JJ+LL 

344 

] 

S(MMfNN)  =  SELMT(I II , J J J ) *SFL*T ( KKK , L LL ) +SELMT  {  M^M, 
L+SELMTI II  I I.JJJJ) 

NNNI 

345 

JJJ=JJJ*1 

346 

LLL=LLL+1 

347 

NNN=NNN«-1 

348 

JJJJ=JJJJ*1 

JJJ=JJJ-4 

LLL=LLL-4 

NNN=NNN-4 
JJJJ  =  JJ.)J-4 

349 

I  I  1=11 1*1 

350 

KKK=KKK+1 

351 

MMM=MMM+1 

352 

II  I  I=IIII  +  1 

353 

GO  TO  371 

354 

00  362  KK=1,4 

355 

MM=I I+KK 

356 

00  360  LL=1,4 

357 

NN=JJ*LL 

358 

S(MM,NN)  =  SELMT( I  II  , J  J  J  )♦  SE  L»H  <  KK  K  ,  L  LL  ) 

359 

JJJ=JJJ«-1 

360 

LLL=LLL*1 
JJJ=JJJ-4 
LLL=LLL-4 

361 

I  I  I  =  I  I  1*1 

362 

KKK=KKK+1 

363 

GO  TO  371 

364 

DO  370  KK=1,4 

365 

MM=I I*KK 

366 

00  369  LL=l,4 

367 

NN=JJ+LL 

368 

S(M^,NN)=SELMT(I II ,JJJ) 

369 

JJJ  =  JJJ«-1 
JJJ=JJJ-4 

370 

111=111+1 

371 

JJ=JJ*4 

372 

I  I  =  I  1+4 

DO  844  I=1,INDEX 
DO  844  J=l,IN0EX 
IF  (OABSISI I ,J) ) .LT.l.D-10)  S(I,l)=0. 

844 

CONTINUE 

USE  80UNOARY  CONDITIONS  TO  MODIFY  S. 
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527  IF    (  (IBNO.GT. 2). OR.UBNO.LT. 2)  )    CO    TO    529 
00    528    1  =  1, M 

II=4*(I-U*N*1 

JJ=ID2 

00    528    KK=1,  INOEX 

S(  II,KK)=0. 

S(KK,U  )=0. 

S( JJ,KK)=0. 

fFK?f  il/llo1  IKK).  AND.  <  Il.EO.D)     SUI,KK)  =  1. 
IF    ( ( JJ.EO.KK) .AND. III. EO. I) )     S(JJ,KK)=1. 

528  CONTINUE 

529  IF  (IBND.LT.3)  GO  TO  531 
DO  530  I=l,M 
II=4*(I-1)*N+1 

JJ=I 1*3 

DO  530  KK=1, INOEX 

00  530  LL=U,JJ 

S(LL,KK)=0. 

S(KK,LL)=0. 

IF  ((LL. EQ.KK)  .AND. (  Il.EO.D  )  S(LL,KK)=1. 

530  CONTINUE 

531  IF  ( (JBND.EQ.2).0R.( JRND.EQ.5) )  GO  TO  543 
GO  TO  533 

543  00  532  1  =  1, M 
II=4*I*N-3 
JJ=ID2 

00  532  KK=1, INDEX 

S< II,KK»=0. 

S(KK,II )=0. 

S( JJ,KK)=0. 

S(KK.JJ)=0. 

IF  UII. EQ.KK). AND. (Il.EO.D)  SUI,KK)  =  1. 

IF  ( (JJ.EO.KK) .AND. ( Il.EO.D )  S(JJ,KK)=1. 

532  CONTINUE 

533  IF  ( ( JBN0.GT.3).0R.( JBN0.LT.3) )  GO  TO  535 
DO  534  1  =  1, M 

II=4*l*N-3 

JJ=II+3 

00  534  KK=1, INOEX 

DO  534  LL=It  ,JJ 

S(LL,KK)=0. 

S(KK,LL)=0. 

IF  (ILL. EQ.KK). AND. < Il.EO.D )  S(LL,KK)=1. 

534  CONTINUE 

535  IF     ( ( JBN0.GT.4).0R.( JBN0.LT.4) )     GO    TO    537 
DO    536    1  =  1, M 

II=4*I*N-2 

JJ=U+2 

DO    536    KK=1, INDEX 

S( II,KK)=0. 

S(KK,II)=0. 

S( JJ,KK)=0. 

S(KK,JJ)=0. 

IF    (  I  II.  EQ.KK)  .AND.  (Il.EO.D  )     SUI,KK)  =  1. 

IF((  JJ.EO.KK).  AND.  (Il.EO.D)     S(JJ,KK)=1. 

536  CONTINUE 

537  IF  (  (KBND.GT. 2). OR.UBND.LT. 2)  )  GO  TO  539 
00  538  1=1, N 

11=4*1-3 

JJ=ID1 

00  538  KK=1, INDEX 

S(  II,KK)  =  0. 

S(KK,II  )=0. 

S( JJ,KK)=0. 

S(KK.JJ)=0. 

IF  ((  II. EQ.KK)  .AND. (  Il.EO.D  )  SUI,KK)  =  1. 

IF  ( (JJ.EO.KK) .AND. ( II. EQ. D )  S(JJ,KK)=1. 

538  CONTINUE 

539  IF  IKBND.LT.3)  GO  TO  541 
DO  540  1  =  1, N 

11=4*1-3 

JJ=II+3 

DO  540  KK=1, INDEX 

DO  540  LL= I  I ,JJ 

S(LL,KK)=0. 

S(KK,LL)=0. 

IF  ( (LL. EQ.KK) .AND. ( Il.EO.D )  S(LL,KK)=1. 

540  CONTINUE 

541  IF  ( (LBND.EQ.2).0R.(LBND.E0.5) )  GO  TO  542 
GO  TO  545 

5*»2  DO  544  1=1, N 

Tl=4*(M-lI*N+4*I-3 

JJ=II+1 

DO  544  KK=1, INDEX 

S( II,KK)=0. 

S(KK,II )=0. 

S(JJ,KK)=0. 

S(KK,JJ)=0. 

IF  ((  II. EQ.KK). AND. (Il.EO.D  )  S(II,KK)=1. 

IF  (  (JJ.EO.KK)  .AND. (  Il.EO.D)  S(JJ,KK(=1. 

544  CONTINUE 

545  IF  ( (LBND.GT. 3). OR.ILBND.LT. 31 )  GO  TO  547 
DO  546  1  =  1, N 

II=4*(M-l)*N+4*I-3 

JJ=I  1*3 

DO  546  KK=1, INDEX 

DO  546  LL=II,JJ 

S<  LL,KK)=0. 

S(KK,LL)=0. 

IF  ( (LL. EQ.KK) .AND. ( Il.EO.D )  S(LL,KK)=1. 

546  CONTINUE 

547  IF  ( (LBND.GT.4).0R.(LBN0.LT.4) )  GO  TO  548 
DP  548  1  =  1, N 

II=4*(M-1)*N*4*I-1 

JJ=II*1 

DO  548  KK=1, INDEX 

S( II,KK)=0. 
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S(KK,II )  =  0. 

SI JJ,KK)=0. 

S(KK.JJ)=0. 

IF  I  I II.EQ.KK) .AND.! Il.EO.l) )  <1CII,KK)  =  1. 

IFII  JJ.EQ.KKKAND. ( Il.EO.l) I  SUJ,KK)=1. 
548  CONTINUE 

IF  ( I1.EQ.2)  GO  Tn  496 
C 

C      SAVE  S  BY  STORING  IN  ARRAY  NAMED  STIFF. 
C 

00  614  1  =  1, INDEX 

00  614  J=1,IN0EX 
614  STIFF! I,J)=S(l,J) 
C 

C       FIND  S  INVERSE  (SINVJ. 
C 

422  DO  426  1=1, INDEX 

423  00  426  J=l, INDEX 

424  SlNVII,J)=0. 

425  IF  I  I.EQ.J)  SINVI I ,J)  =  1 

426  CONTINUE 

427  DO  443  1=1, INDEX 

428  K=I 

429  IF  (I-INDEX)  430,437,430 

430  IF  (S(I.I)-FPS)  431,432,437 

431  IF  (-SUtl)-EPS)  432,432,437 

A  1  *>   |(  —  «  +  1 

433  DO  435  J=l. INDEX 

4  34  SINVI I,J)=SINV{ I ,J)+SINVIK,JI 

435  S(  I, J)  =  S(  I  ,J)+S(K,J) 

436  GO  TO  430 

437  DIV=S(  1,11 

438  00  440  J=l, INDEX 

439  SINVI I,J)=SINV(l ,J)/OI V 

440  S( I,J)  =  S( I  ,J)/OI V 

441  00  448  L  =  l, INDEX 

442  DELT=S(L,I I 

443  IF  (DABSIDEl Tl-FPS)  448,443,444 

444  IF  (L-I)  445,443,445 

445  00  447  J=l,  INDEX 

446  SINVIL,J)=SINV<L,J)-SINVI I ,J)*OELT 

447  S(L, J)=S(L,J)-S( I ,J)*DELT 

448  CONTINUE 

00  845  1=1, INDEX 
DO  845  J=l,  INDEX 

IF  (DARSI SINVI  I, J)  I. IT. 1.0-10)  SINVII,J)=0. 
845  CONTINUE 

C  THIS  PORTION  OF  THE  PROGRAM  PRODUCES  THE  RECTANGULAR  "LATF  ELEMFNT 

C  STABILITY  MATRIX  RY  A  PROCEDURF  SIMILAR  TO  THAT  FOR  THP  FLEMFNT 

C  STIFFNESS  MATRIX 

C 

C  ESTABLISH  LOAD  CONDITION  MATRIX  C. 

C 

449  C(  1,  1)=C2 

450  CI  1,2)=C3 

451  C(2,1)=C3 

452  CI2,2)=C1 
C 

C  ESTABLISH  MATRIX  H 
C 

453  00  459  1=1,16 

454  HCOEFI 1,1 )=PCLY2( I ) 

455  HC0EFI2.I )=PCLY1I I ) 

456  HXEXPI 1,1 )=P0LY1 II ) 

457  HXEXP(2,I )=PCLY3I I ) 

458  HYEXPI 1, I )=PCLY4( I ) 

459  HYEXPI 2,1 )=PCLY?< I ) 
C 

C  THE  FOLLOWING  PART  OF  THF  PROGRAM  COMPUTES  INTEGRAL  OVFR  X  FROM  0 

C  TC  EX,  INTEGRAL  OVER  Y  FROM  0  TO  EY,  OF  H  TRANSPOSE  *  c  *  H. 

C  THE  RESULT  IS  NAMFD  HINT. 
C 

460  DO  463  1=1,16 

461  DO  468  J=l  ,16 

462  SUM1=0. 

463  DO  467  K=l,2 

464  SUM2=0. 

465  00  466  L=l ,2 

4  66  SUM2=SUM2*HCOEF(K,  I  )  *C  (  K  ,  L  )  *HC  OEF  (  L  ,  J  )  *  (  l./I(HXFX°<K,I)-»-HXFXP(L,J) 
l  +  l.  1*1 HYEXPI K,  I 1+HYFXPIL, J>  +  1 .  )  )  )*( FY *♦ ( HYF X  PI K , I  UHYEXPII  , J  )  ♦  1  .  )  ) 
2*IEX**(HXEXP(K,I  )*HXrXP(L,J)  +  l.)) 

467  SUM1=SUM1*SUM2 

468  HINTI I,J)=SUM1 
C 

C  FIND  ELEMENT  STABILITY  MATRIX  (8EL"T).  BFLMT=AINV  TRANSPOSE  * 

C  HINT  *  AINV. 

C 

476  DO  484  1=1,16 

477  DO  484  J=l,16 

478  SUM1=0. 

479  DO  483  K=l,16 

480  SUM2=0. 

481  DO  4  82  L  =  1  ,16 

482  SUM2=SUM2*HINTIK,L)*AINV(L,J) 

483  SUM1=SUM1+AINVIK,I )*SUM2 

484  RELMTI  I  ,  J  >=SUMl 
C 

C  GENERATE  ASSEMBLED  STABILITY  MATRIX  AND  APPLY  BOUNDARY  CONDITIONS 

C  RY  USING  SAME  PORTION  OF  PROGRAM  USED  Fnp  STIFFNESS  MATRIX. 

C 

494  11=2 

495  GO  TO  209 
4"6  CONTINUE 

C 

C      THIS  PORTION  FINDS  THE  RUCKLING  COFFFICIENT  AND  DISPLACFMFMT 

C      VECTOR  BY  MATRIX  ITERATION. 

C 
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L=l 

MM=1 
C 

C      ESTABLISH  TRIAL  VECTOR  OF  ELEMENTS  RETWEEN  -1  AND  *1  TAKFN  FROM 
C      TABLE  OF  RANDOM  DIGITS 

11  =  1 
JJ=10 

497  DO  739  1=1,10 

READ  (5,740)  (W( JJ  ,J=I I,JJ) 
11=11*10 

741  FORMAT  tlHO, 19HTRI AL  VECTOR,  F IRST , 14, 1 X , 17HELEMENTS  OF  ARRAY) 
11  =  1 

JJ=10 

DO  498  1=1,10 

WRITE  (6,742)  ( W( J ) , J=I I , J J) 

1 1=1  1*10 

JJ=JJ*10 

498  CONTINUE 

742  FORMAT  (10F6.2) 
C 

C      MULTIPLY  W  TRANSPOSE  *  S  *  W  =  WTSW. 
C 

700  WTSW=0. 

DO  701  1=1, INDEX 
00  701  J=l, INDEX 

701  WTSW=WTSW*W<  I)*STIFF( I,J)*W( J) 
C 

C      NORMALIZE  W. 
C 

DO  702  1=1, INDEX 

702  W<  I)=W( I)/DSQRT(WTSW) 
IF  (MM.EO.l)  GO  TO  713 

C 

C      AFTER  FINAL  I TER AT  I  ON, WR I TE  TRANSVERSE  DISPLACEMENTS  IN  MAP  FORMAT 

C 

722  IF  (C3.EQ.0.)  GO  TO  721 
DO  716  1=1, INDEX 

IF  (MM.E0.2)  W(I )=WSAVE( I)+W(t ) 
IF  (MM.E0.3)  W(l )=2.*WSAVE(I)-W(I ) 
716  CONTINUE 
721  WRITE  (6.608) 

608  FORMAT  ( 1H1, 22X, 24HTRANSVERSE  DISPLACEMENTS) 
WRITE  (6.609) 

609  FORMAT  < 1H0, 24X, 20HY  INCREASES  TO  RIGHT/) 
WRITE  (6,610) 

610  FORMAT  (21X.26HX  INCREASES  TOWARDS  BOTTOM/) 
K=l 

L=4*N-3 

00  612  J=1.M 

WRITE  (6,611)  (W(I ),l=K,L,4) 

611  FORMAT  (8(D11.4,5X)I 
WRITE  (6,613) 

613  FORMAT  ( 1H  //////////) 
K=L+4 

612  L=L*4*N 

IF  (C3.E0.0.)  GO  TO  712 

MM=MM+ 1 

IF  (MM.E0.3)  GO  TO  722 

IF  (MM.E0.4)  GO  TO  712 
C 

C      MULTIPLY  SINV  *  8  *  W( NORMAL  I  ZED )  =  WNEXT 
C 

713  DO  705  1=1, INDEX 

SUM2=0. 

DO  704  J=l, INDEX 

SUM1-0. 

DO  703  K=l, INDEX 

703  SUM1=SUM1+B( J,K)*W(K) 

704  SUM2=SUM2*SINV(I ,J)*SUM1 

705  WNEXT( I )=SUM2 
C 

C      MULTIPLY  WNEXT  TRANSPOSE  *  S  *  WNEXT  =  EIGEN 
C 

EIGEN=0. 

DO  706  1=1, INDEX 

DO  706  J=l, INDEX 

706  EIGEN=EIGEN+WNEXT( I ) *STIFF ( I , J )*WNEXT ( J) 

C      FIND  BUCKLING  COEFFICIENT. 
C 

AK=1./(DSQRT(E!GEN)*(3.14159265*AR)**2) 
WRITE  (6,606)  L.AK 
606  FORMAT  < 1H0, 9H I TER AT  I  ON, I  3 ,4H ,  K  =  ,F9.5) 

C      TEST  FOR  20  ITERATIONS. 
C 

IF  (L.E0.20)  GO  TO  710 
C 

C      TEST  FOR  CONVERGENCE. 
C 

IF  (L.EO.l)  SAVF=AK 

IF  (L.EO.l)  GO  TO  708 

IF  (DABS(OABS(AK)-DABS(SAVE))/DABS( AK).LT.. 00001)  GO  TO  711 

SAVE=AK 
708  L*L*l 
C 

C      SAVE  W  (NORMALIZED) 
C 

00  715  1=1, INDEX 
715  WSAVE( I )=W(I  ) 
C 
C      USE  RESULTING  VECTOR  AS  NEXT  TRIAL  VECTOR. 


C 


714  DO  709  1  =  1, INDEX 
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709  W(  I  )=WNEXT(  I  ) 
GO  Tn  700 

710  WRITE  (6,522) 

522  FORMAT  (1H0,VM-I20  ITERATIONS  WITHOUT  CONVFR  Of  no  F  ) 
GO  TO  712 

711  WRITF     (6.525) 

525     FORMAT     <  lH0,48HtUJ0KI   ING    PARAMrTFR    HAS    CONVEROEO    TO     .DTI     Profr-.iT) 
MM=2 
GO    TO    722 

712  CONTINUE 
FNO 


ASPtOT  RATIU=  1.00000    P01SSMN  PATIO  =  0.0     NO.  ELEMENTS  IN  X  =  4    NO.  ELEMENTS  IV  Y  =  4 
3UUNDAPY  CONDITIONS,  1=FKCF,  2=PINMEU,  3=CLAMPE0,  4=SYM.,  5  =  ANTI-<rYM. 

Y=0,flJGF  CONDITION  3    Y=B,EnGE  CONDITION  4    X  =  0,FD';E  CONDITION  2    X=AfFDGE  CONDITION  c 
^FLATIVE  LOAD  S17ES   NX=  1.00   NY=  0.0    NXY=  0.0 

TgTAL  VECIORf  FIRST  10Q  ELEMENTS  fT  ^^SY 

0.07       0.52  -0.72       0.17  -0.13       0.2^    -0.96  0.68    -0.54  0.6* 

-0.09    -0.33  -0.65       0.72  0.8w       0.P1    -0.17  -0.54       0.22  0.75 

-O.CJ    -0.11  0.60    -0.79  0.71    -0.C9       0.63  0.1."    -0.17  0.5^ 

0.22       0.14  0.91       0.14  0.13       0.58       0.50  -0.93       0.60  -0.63 

U.16    -0.47  0.27    -0.7*  -0.55        0.«5        0.17  0.98    -0.25  0.04 

0.3'^       0.42  -0.98       0.97  -0.60    -0.04    -0.18  0.30       0.15  -0.61 

0.31    -0.65  -0.57    -0.7*  -0.76       0.97       0.21  0.98    -0.21  0.21 

0.73    -C.33  0.66    -0.68  -0.19    -0.4°    -0.76  0.55    -0.33  -0.1'7 

0.*3       0.34  0.99       0.95  0.38    -0.?Q       0.41  -0.95       0.98  0.75 

-O.dl     -0.74  0.34       0.9Q  0.55    -0.31        0.94  -0.06    -0.69  -0.23 


I TERATION 

1, 

K  = 

104. 

51617 

I TERATION 

2, 

K  = 

26. 

,53432 

ITFRATICN 

3, 

K  = 

13. 

,08651 

I TERAT  ION 

4, 

K  = 

8. 

81542 

ITERATION 

5, 

K  = 

7. 

,890  34 

ITERATION 

6, 

K  = 

7. 

,72830 

I TERATION 

7, 

K  = 

7. 

70094 

ITERATION 

8, 

K  = 

7. 

69634 

ITFSATICN 

9, 

K  = 

7. 

6955S 

I TERATION 

10, 

K  = 

7. 

£.9*43 

ITFRATION 

U. 

K  = 

7. 

69541 

■SUCKLING    PARAMETER    HAS    CONVERGED    TO    .001     PERCENT 
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TRANSVERSE  DISPLACEMENTS 
y  INCREASES  TO  RIGHT 
X  INCREASES  TOWARDS  BOTTOM 
0.0  0.0  0.0  0.0  0.0 


0.0  -0.1212D-01      -0.3320D-01      -0.50040-01      -0.56320-01 


0.0  -0. 17150-01      -0.46970-01      -0.7081D-01      -0.79690-01 


0.0  -0.1214D-01      -0.33230-01      -0.50090-01      -0.56380-01 


0.0  0.0  0.0  0.0  0.0 
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